module suplmod
implicit none
contains
subroutine suplmf(m,T,n,y,theta,lmk,lms,suplm,kt)
use linear_operators
use matrixfuns
use datahadler
integer, intent(in):: m,T,n
double precision, intent(in):: theta(m),y(T,n)
double precision, intent(out):: lmk,lms,suplm,kt
integer i,j
double precision, parameter:: c9=.999d0
double precision, parameter:: pi=3.14159265358979d0
double precision lamb(T),delta(n),phi(n),c(n),phi0(n),phi1(n),phi2(n),&
                &alpha1,alpha2,llik,sigmat(n,n),gamma(n),gammat(n,n),&
				&cc(n,n),epst(1:T,1:n),isigmat(n,n),igammat(n,n),alpha0&
				&,fttm,wttm,lmki,lmsi(n),asigt,isigmatm(n,n)

delta=theta(1:n)
c(1)=1.0d0
c(2:n)=theta(n+1:2*n-1)
alpha0=theta(2*n)
phi0=vassign(n,theta(2*n+1))
alpha1=c9*(dsin(theta(2*n+2))**2.0d0)
alpha2=(c9-alpha1)*(dsin(theta(2*n+3))**2.0d0)

forall(i=1:n)
phi1(i)=c9*(dsin(theta(2*n+4))**2.0d0)
end forall
forall(i=1:n)
phi2(i)=(c9-phi1(i))*(dsin(theta(2*n+5))**2.0d0)
end forall

forall(i=1:n,j=1:n)
	cc(i,j)=c(i)*c(j)
end forall

forall(i=1:T,j=1:n)
	epst(i,j)=y(i,j)-delta(j)
end forall


gamma=phi0/(1.0d0-phi1-phi2)


gammat=diag(gamma)

lamb(1)=alpha0/(1.0d0-alpha1-alpha2)

igammat=diag(1.0d0/gamma)
sigmat=(cc*lamb(1))+gammat

isigmat=igammat-lamb(1)*(((igammat.x.cc).x.igammat)&
     &/(lamb(1)*dot_product(c,matmul(igammat,c))+1.0d0))

isigmatm=isigmat

asigt=dot_product(epst(1,:),matmul(.i.sigmat,epst(1,:)))
lmki=.25d0*(asigt**2.0)-.5d0*(n+2.0d0)*asigt+.25d0*n*(n+2.0d0)
lmsi=(asigt-n-2.0d0)*epst(1,:)


do i=2,T,1
	wttm=lamb(i-1)/(1.0d0+lamb(i-1)*dot_product(c,c/gamma))
	fttm=wttm*dot_product((c/gamma),epst(i-1,:))
	lamb(i)=alpha0&
	       &+alpha1*(((fttm)**2.0d0)+wttm)+alpha2*lamb(i-1)


	gamma=phi0+phi1*(((epst(i-1,:)-c*fttm)**2.0d0)+(c*c)*wttm)+phi2*gamma

	gammat=diag(gamma)

	igammat=diag(1.0d0/gamma)

	sigmat=(cc*lamb(i))+gammat
	isigmat=igammat-lamb(i)*(((igammat.x.cc).x.igammat)&
	     &/(lamb(i)*dot_product(c,matmul(igammat,c))+1.0d0))

	asigt=dot_product(epst(i,:),matmul(isigmat,epst(i,:)))
	lmki=lmki+.25d0*(asigt**2.0)-.5d0*(n+2.0d0)*asigt+.25d0*n*(n+2.0d0)
	lmsi=lmsi+(asigt-n-2.0d0)*epst(i,:)

end do
lmk=(2.0d0/(n*(n+2.0d0)))*lmki*lmki/T
lms=(1.0d0/(2.0d0*(n+2.0d0)))*dot_product(lmsi,matmul(isigmatm,lmsi))/T
suplm=lmk+lms
if (lmki>0) then
	kt=suplm
else
	kt=lms
end if

end subroutine suplmf

subroutine givebin(m,T,n,y,theta,bin)
use linear_operators
use matrixfuns
use datahadler
integer, intent(in):: m,T,n
double precision, intent(in):: theta(m),y(T,n)
double precision, intent(out):: bin(n)
integer i,j
double precision, parameter:: c9=.999d0
double precision, parameter:: pi=3.14159265358979d0
double precision lamb(T),delta(n),phi(n),c(n),phi0(n),phi1(n),phi2(n),&
                &alpha1,alpha2,llik,sigmat(n,n),gamma(n),gammat(n,n),&
				&cc(n,n),epst(1:T,1:n),isigmat(n,n),igammat(n,n),alpha0&
				&,fttm,wttm,lmki,lmsi(n),asigt,isigmatm(n,n),mst(n),mkt&
				&,sigmabar(n,n)

delta=theta(1:n)
c(1)=1.0d0
c(2:n)=theta(n+1:2*n-1)
alpha0=theta(2*n)
phi0=vassign(n,theta(2*n+1))
alpha1=c9*(dsin(theta(2*n+2))**2.0d0)
alpha2=(c9-alpha1)*(dsin(theta(2*n+3))**2.0d0)

forall(i=1:n)
phi1(i)=c9*(dsin(theta(2*n+4))**2.0d0)
end forall
forall(i=1:n)
phi2(i)=(c9-phi1(i))*(dsin(theta(2*n+5))**2.0d0)
end forall

forall(i=1:n,j=1:n)
	cc(i,j)=c(i)*c(j)
end forall

forall(i=1:T,j=1:n)
	epst(i,j)=y(i,j)-delta(j)
end forall


gamma=phi0/(1.0d0-phi1-phi2)


gammat=diag(gamma)

lamb(1)=alpha0/(1.0d0-alpha1-alpha2)

igammat=diag(1.0d0/gamma)
sigmat=(cc*lamb(1))+gammat

isigmat=igammat-lamb(1)*(((igammat.x.cc).x.igammat)&
     &/(lamb(1)*dot_product(c,matmul(igammat,c))+1.0d0))

isigmatm=isigmat

asigt=dot_product(epst(1,:),matmul(.i.sigmat,epst(1,:)))
lmki=.25d0*(asigt**2.0)-.5d0*(n+2.0d0)*asigt+.25d0*n*(n+2.0d0)
lmsi=(asigt-n-2.0d0)*epst(1,:)
sigmabar=sigmat

do i=2,T,1
	wttm=lamb(i-1)/(1.0d0+lamb(i-1)*dot_product(c,c/gamma))
	fttm=wttm*dot_product((c/gamma),epst(i-1,:))
	lamb(i)=alpha0&
	       &+alpha1*(((fttm)**2.0d0)+wttm)+alpha2*lamb(i-1)


	gamma=phi0+phi1*(((epst(i-1,:)-c*fttm)**2.0d0)+(c*c)*wttm)+phi2*gamma

	gammat=diag(gamma)

	igammat=diag(1.0d0/gamma)

	sigmat=(cc*lamb(i))+gammat
	isigmat=igammat-lamb(i)*(((igammat.x.cc).x.igammat)&
	     &/(lamb(i)*dot_product(c,matmul(igammat,c))+1.0d0))

	asigt=dot_product(epst(i,:),matmul(isigmat,epst(i,:)))
	lmki=lmki+.25d0*(asigt**2.0)-.5d0*(n+2.0d0)*asigt+.25d0*n*(n+2.0d0)
	lmsi=lmsi+(asigt-n-2.0d0)*epst(i,:)

	sigmabar=sigmabar+sigmat
end do
sigmabar=sigmabar/T
mkt=lmki
mst=lmsi
bin=(n/(4.0d0*mkt))*matmul(.i.sigmabar,mst)
end subroutine givebin

end module suplmod